Black Ash Wetland Response to Future Climate Conditions

Abstract Future climate scenarios will likely lead to reduced connectivity between black ash headwater wetlands and stream networks.

Black Ash Wetland Response to Future Climate Conditions

Introduction

Shifts in temperature, and precipitation (timing, frequency, and quantity) associated with climate change are already impacting global wetland ecology and will intensify into the future (Burkett and Kusler 2000; Moomaw et al. 2018). In general, the Great Lakes region of North America (Figure 1) is expect to see declining summer precipitation while total annual precipitation stays stable or increases (Byun and Hamlet 2018; Hayhoe et al. 2010). These projected changes result in increased precipitation in the spring and winter months. Winter precipitation will also experience a phase change, seeing a reduction in the proportion of precipitation as snow and increased occurrence of rain-on-snow melt events in the winter and spring (Notaro, Bennington, and Vavrus 2015). Locally conditions will vary and our study area (see below) is projected in numerous downscaling scenarios across multiple GCMs to see a smaller than regional decline or even slight increases in summer precipitation (LOCA sitation, https://climatetoolbox.org/tool/Climate-Mapper, Byun and Hamlet 2018, fig. 7). Taken together the shift in precipitation timing, the reduction of total snowmelt, and an earlier snowmelt will reshape the annual hydrologic budget. These changes will alter wetland hydroperiods, with water availability increasing in the winter and spring and decreasing in the summer and early fall. The impact of reduced or steady precipitation in the summer will be compounded by increased evaporative demand (Byun and Hamlet 2018). Summer temperatures are expected to increase 2-8C° (Byun and Hamlet 2018; Hayhoe et al. 2010), and without a commenserate increase in summer precipitation will lead to more frequent drought conditions and increased water stress on ecosystems.

Black ash (Fraxinus nigra Marsh) is an important hardwood component of many forested wetlands in the northern United States and southern Canada. Wetlands with a large or majority black ash component face the same climate change impacts as other wetlands in the region. They also face loss of the major canopy species due to an invasive insect, Emerald Ash Borer (Agrilus planipennis Fairmaire, EAB). Emerald ash borer was first found in the United States in southeast Michigan in 2002 (Haack et al. 2002). It is known to infest and cause high mortality in all ash native to North America (Herms and McCullough 2014). As of this writing it is present in 35 US states and 5 Canadian provinces. High mortality and the importance of ash in regional forested wetlands has led to research on the impacts of their loss and strategies to mitigate those impacts. Iverson et al. (2016) evaluated the potential replacement species for black ash in the context of habitat availability, species migration, and replacement species susceptibility to climate change impact, providing a useful resource for climate-informed species replacements for black ash. However, black ash grows in a range of geomorphic settings and local site conditions can have a strong influence on hydrology, forest structure, and plant community among wetlands (Kolka et al. 2018). Considering both Iverson et al. (2016) and Kolka et al. (2018) we can identify four necessary components to evaluate mitigation efforts in black ash wetlands: 1) climate-informed species selection, 2) species tolerance of local site conditions, 3) site conditions following EAB infestation, and 4) site conditions in a future climate.

Previous and ongoing work in Michigan and Minnesota has assessed mitigation in light of three of the four components of the combined impact of climate change and EAB on black ash wetlands. Researchers in Michigan (Bolton et al. 2018), Minnesota (Looney et al. 2015), and Wisconsin (Bolton et al. 2018), planted seedlings to evaluate potential canopy species at the wetland rather than landscape level. The plantings in Bolton et al. (2018) and Looney et al. (2015) took place under simulated EAB infestation, where the seedlings were subjected to adverse conditions due to the hydrologic impact of the loss of black ash (Slesak et al. 2014; Van Grinsven et al. 2017), as well as the increased competition from herbaceous growth under increased light conditions (Looney et al. 2016, 2017; Davis et al. 2017).

The researchers addressed 3 of the 4 components necessary for mititgation efforts in black ash wetlands: 1) some planted species were at the northern edge of their current range and evaluated 2) under present-day site conditions and 3) under simulated-EAB site conditions. What is not possible in field studies is addressing the interaction of EAB and climate change on site conditions. Just as we have seen EAB impacts cascade through black ash ecosystems affecting hydrology, plant communities, and nutrient cycling, we can expect the future climate-driven changes to hydrology to result in similar cascades. Focusing on the future hydrologic characteristics of these wetlands will inform us on the impacts to other functions, as hydrology is the critical control on wetland ecosystems (Brinson 1993). In order to understand future hydrologic conditions we have derive wetland water level models and evaluate under future climate scenarios.

Simulated daily time-step weather conditions are available for various future climate scenarios to serve as model inputs, providing realizations of 30-year climate periods. However, previous studies have demonstrated high-interannual variation in seasonal wetland water level drawdown and rebound in black ash wetlands (Van Grinsven et al. 2017). Estimating the behavior of a highly variable system with 30 years of daily climate projects could lead to biased or high-variance results making it difficult to draw proper conclusions about the combined impacts of climate change and EAB. A 30-year period provides a representation of normal climate variation, but the weather sequences underlying a certain climate can take many more forms than the observed (or simulated) 30 annual weather sequences. Therefore, simulation studies of alternative scenarios or future conditions require more weather sequences to quantify the most likely response and the distribution of possible responses. Stochastic weather generators (SWGs) provide a tool for generating synthetic time series of weather that simulate conditions under an observed or projected climatology (D. S. Wilks and Wilby 1999). SWGs have also been used as a downscaling technique for global or regional climate models (Daniel S. Wilks 2012; Verdin et al. 2015, 2019). One class of SWG known as Richardson (Daniel S. Wilks 2012) are built from parametric representations fit to climate conditions using multiple generalized linear models and varied statistical distributions. This approach can be implemented using any observed (or simulated) weather series that define particular climatic conditions. This characteristic makes them well suited to compare multiple climate scenarios without the overhead cost of executing complete GCM or regional climate model (RCM) runs.

We have performed simulation experiments combining observed wetland hydrology and synthetic weather sequences to quantify the interactions of EAB and future climate scenarios. We developed wetland hydrology models for current canopy conditions (ash-dominated forested wetlands), non-forested conditions (herbaceous and shrub/scrub wetlands), and potential future canopy conditions (forested wetlands under current co-dominant species). We evaluated each class of model under two potential future climate scenarios for the end of the 21st century (2070-2099). The two future scenarios are defined by a less sensitive GCM (generally projects smaller magnitude regional climate impacts) under the moderate relative concentration pathway (RCP 4.5), and a more sensitive GCM (projects larger magnitude regional climate impacts) under the relative concentration pathway (RCP 8.5) (van Vuuren et al. 2011). The two-scenario ‘bookend’ approach provides a range of potential future conditions as opposed to a multi-model ensemble approach which masks some of the uncertainty in potential future conditions (Swanston et al. 2016). For an unbiased comparison between current and future conditions both GCMs are used to evaluate the wetland models under historic (1980-2009) climate conditions.

We expect that the interaction of hydrologic impacts of EAB and climate change will result in a tempering of the two individual impacts in this region. While simulation of post-EAB conditions have lead to increased water levels and reduced drawdown rates in growing season, future climate conditions in the region will result in reduced water availability during that same period. These two opposing drivers should result in some moderation to both impacts.

Methods

Study Area & Data

Wetland water levels measured at fourteen black ash-dominated wetlands in the western Upper Peninsula of Michigan, USA were used to develop and evaluate our wetland hydrologic models (Fig. 1). The wetlands range in size from 0.29 to 1.19 ha and 30–80% of the basal area consists of black ash with histosol soils over a confining layer located at an average depth of 118.8 cm (Davis et al. 2017; Van Grinsven et al. 2017). The region has average minimum and maximum annual temperatures of -11.3 °C and 18.2 °C and an average annual precipitation of 101 cm for the climate period of 1980-2009 at the Bergland Dam (46°35’13”N, 89°32’51”W) meteorological station (Arguez et al. 2012). Wetland water levels were continuously monitored in 2” inner-diameter driven wells from 2012-2020 and logged every fifteen minutes using Solinst Levellogger Junior pressure transducers (Solinst, Ontario, CA), with more details available in (Van Grinsven et al. 2017). Barometric compensation was performed using data from Solinst Levellogger Junior pressure transducers deployed at a subset of study wetlands. Compensated water levels were corrected for temperature differentials as in Shannon, et al (In Review). Following an initial control period of two growing seasons, two-thirds of the wetlands in the study were treated to simulate the impacts of an EAB infestation. At one-third of the sites all ash trees greater than 1” in diameter were girdled and at the other half of the treatment sites all ash trees greater than 1” in diameter were hand-felled and left on site.

Daily precipitation and daily minimum/maximum temperatures from existing meteorological stations were used as inputs to the wetland water level models described below. From these input drivers we derived daily solar radiation, PET, precipitation as snowfall and rain, and snowmelt. Precipitation records were retrieved from the National Centers for Environmental Information Hourly Precipitation Dataset (“Hourly Precipitation Data (HPD) Network, Version 2.R2” 2021) using the stations USC00201088 (Bruce Crossing, MI), USC00204328 (Kenton, MI), USC00475352 (Mercer Ranger Station, WI), USC00206215 (Ontonagon, MI), USC00476398 (Park Falls, WI), USC00476518 (Phelps, WI), USC00476939 (Rainbow Reservoir Lake, Tomahawk, WI), USC00477140 (Rice Reservoir, Tomahawk, WI), and USC00208680 (Watersmeet Fish Hatchery) and summed to daily values. Daily minimum and maximum temperatures were taken from the Global Historical Climatology Network (GHCND) dataset (Menne, Durre, Korzeniewski, et al. 2012; Menne, Durre, Vose, et al. 2012) for the same stations, which were collocated stations for HPD and GHCND. Where data were missing, values were filled using inverse distance weighting with data retrieved from the Mesowest (Horel et al. 2002) network stations BPLM4 (Baraga Plains, MI), KTNM4 (Kenton, MI), PIEM4 (Pelkie, MI), WKFM4 (Wakefield, MI), WMTM4 (Watersmeet, MI). Solar radiation data were also retrieved from the listed Mesowest stations. The solar radiation and temperature data were used to fit the Bristow-Campbell method coefficients to estimate solar radiation from latitude and daily temperature range using the PIEM4 Mesowest site (Bristow and Campbell 1984; Bojanowski 2016). Bristow-Campbell solar radiation was used to calculate potential evapotranspiration (PET) via the modified Hargreaves-Samani equation (Hargreaves and Allen 2003). Precipitation was partitioned into snowfall, rain, and snowmelt inputs using the CemaNeige snow accounting routine (SAR) (Valéry, Andréassian, and Perrin 2014). The CemaNeige SAR is temperature-index based and accounts for accumulation, snowpack cold content, and snowmelt through a thermal state weighting coefficient and a degree-day melt coefficient. The CemaNeige SAR is implemented in the R package airGR (Coron et al. 2017).

Figure 1: Map of study site and meteorological stations used for data retreival. Coordinates in meters, UTM Zone 16N.

Wetland Hydrology Models

The link between hydrology drivers (rainfall, snowmelt, PET) and daily water level response has been shown to be non-linear and varying with stage (Loheide II 2008; White 1932). In wetland systems the relationship between the driver magnitude and response magnitude has been termed the ecosystem specific yield (E_Sy) (McLaughlin and Cohen 2014) (see Section #### in Chapter #### on Pressure Transducer Uncertainty for more information about the development of E_Sy over time). E_Sy has previously been empirically derived as the ratio between precipitation inputs and water level rise (\(\frac{P}{\Delta\;WL}\)) (McLaughlin and Cohen 2014). The relationship between empirical E_Sy and water level can then be modeled to provide a continuous estimate of E_Sy. Models fit to the E_Sy~Water Level relationship include exponential (McLaughlin and Cohen 2014; Watras et al. 2017), quadratic (McLaughlin and Cohen 2014), and step-wise regression (McLaughlin and Cohen 2014). Identifying E_Sy using the rainfall-rise ratio can be unworkable in the presence of confounding hydrologic variation such as surface water connectivity and low-frequency seasonal changes in water availability (Watras et al. 2017; Zhu et al. 2011). Both of these factors were present in our study wetlands and no clear relationship existed to model E_Sy with the above model forms, requiring an alternative approach.

As an alternative we fit an inverse analog to \(\frac{P}{\Delta\;WL}\), deriving E_Sy from the ratio of cumulative water availability and water level. We began by fitting a quadratic curve to the relationship between a year-to-date water availability index from the beginning of the growing season to the point of minimum wetland water level (Figure 2), where

\[WA_{YTD}\;=\;P_{YTD} + M_{YTD} - PET_{YTD},\]

\(WA_{YTD}\), \(P_{YTD}\), \(M_{YTD}\), and \(PET_{YTD}\) are year-to-date water availability, rainfall, snow melt, and potential evapotranspiration, respectively. Empirical E_Sy was then considered as the first derivative of the quadratic curve, which has the form \(\frac{\Delta\;WL}{\Delta\;WA}\). This approach resulted in an asymptotic relationship between E_Sy and water level, suggesting agreement with the exponential forms in McLaughlin and Cohen (2014) and Watras et al. (2017). Defining E_Sy as \(\frac{\Delta\;WL}{\Delta\;WA}\) provided additional information to fit the relationship between E_Sy and water level as we were not limited to days with rainfall. Models for E_Sy were fit using the year with the greatest water level drawdown for each wetland to capture the widest range of E_Sy variation. Each wetland had the same basic form of E_Sy function. The functions were fit using a single hierarchical model with each of the coefficients allowed to vary independently within sites. Models were fit using the brms package in R (Bürkner 2017, 2018). The asymptotic structure of the E_Sy function requires that a lower bound be placed on E_Sy predictions to avoid values less than or equal to zero. The minimum value of E_Sy was allowed to vary by wetland and was fit along with other model parameters as described below.

Figure 2: Derivation of ecosystem specific yield for a single study wetland. Panel A shows the quadratic relationship between wetland water level and water balance (total liquid inputs minus potential evapotranspiration) for the drawdown period. Panel B shows the relationship between the derivative of the fitted line from Panel A, representing an emprical estimate of E_Sy and wetland water level.

Wetland water levels were simulated on a continuous basis using daily inputs. For each daily step, wetland water level change was determined as a function of the current water level, daily rainfall (R), daily PET, daily snowmelt (M), and estimated streamflow (Q) (Equation EQAAA; Table 1). Maximum wetland water levels were estimated from the the control period as the mode of the wetland water level record. Streamflow was assumed to occur whenever wetland water levels were at or above the maximum water level. This is similar to the approach in McLaughlin et al. (2019) where wetland surface water connectivity was determined from wetland water level records. The current water level was used to estimate \(\hat_{E}~Sy~\) which was then used as a multiplier for rainfall, PET, and snowmelt components. In addition each driver (R, PET, M, Q) had a coefficient (fitting described below). Snowmelt and precipitation were also fit with a first order autoregressive filter to simulate slow flow contributions to wetland water levels. Within a single day only the larger of PET or P was used as driver of water level change, accounting for suppressed transpiration in from wet leaves.

\[ \begin{align*} & WL_{t+1} = WL_{t} + f_{E_{Sy}}(WL_{t-1})*(\hat{\beta}_DD + \hat{\beta}_MM) + \hat{\beta}_Q(WL_{t} - WL_{max})\\ & \hat{\beta}_DD\;=\left\{R \le PET;\;\;\hat{\beta}_{PET}PET\atop{R > PET; \;\; \hat{\beta}_RR}\;\;\;\;\;\;\;\;\;\;\right.\\ & f_{E_{Sy}}(WL_{t-1})\;=\;a - (a - b) * e^{(c * WL_{t-1})} \\ & M_t = M_t + \phi_M\; M_{t-1} \\ & R_t = R_t + \phi_R\;R_{t-1} \end{align*} \]

Training data for wetland water level models was selected for each wetland as the year with the greatest water level drawdown within the control period. The remaining control-period years were used for wetland model evaluation. Wetland parameters were estimated using the L-BFGS-B algorithm via the R optim() function (R Core Team 2019) minimizing the weighted root-mean square error of the predictions. Weights were applied asymmetrically where days with observed water levels greater than or equal to WL_max were asszigned the original equal weight (1/n). As observed water levels decreased below WL_max weights increased as the square of the difference between observed water level and WL_max. To model impacted wetlands that have not recovered to a fully forested state the observed treatment data from the same site were used to reparameterize the control condition model. Treatment-period training data were also selected to maximize water level drawdown in the training set. The wetland models were reparameterized for only \(\hat{\beta}_{PET}\), \(\hat{\beta}_{R}\) terms. Finally we simulated reforested black ash wetlands under one set of potential future forest composition. For the future forest conditions we assumed a mix of the current co-dominant species, red maple and yellow birch, would become established with similar stand basal area. Future forested simulations were performed using the same parameters as the control period with a reduction of \(\hat{\beta}_{PET}\). The reduction in \(\hat{\beta}_{PET}\) was a function of water level and current proportion of site basal area as black ash:

\[\hat{\beta}_{PET-Adj} = \hat{\beta}_{PET}*[(1 - BA_{ash}) + BA_{ash} * [1.45077 - 0.05869 * (WL-WL_{max})]^{-1}]\] This equation is derived from previous work showing a water level-depended difference in sap flux between black ash and current co-dominant species (Shannon et al. 2018).

Because individual models were trained for each wetland there are no training data available for non-forested conditions under the field-study control sites. Therefore we fit wetland water level models to only of 8 of the field-study wetlands: the 8 treatment (girdle and ash cut) wetlands. Comparisons of the combined and separate impacts of EAB and climate change were performed against modeled baselines rather than field-study observed conditions. Limiting the number of sites and comparing to modeled future simulated baselines with each alternative vegetation condition ensures that comparisons are not biased by the number of hydrology of sites within each group. Control conditions are represented by the modeled baselines, avoiding the potential of identifying model artifacts as significant when comparing observed and modeled results.

Future Climate Conditions

Wetland models were run under simulated current (1980-2009) and future (2070-2099) conditions. Future conditions were evaluated under the RCP 4.5 and RCP 8.5 representative concentration pathway scenarios (van Vuuren et al. 2011). Daily climate projections for all scenario-period combinations were taken from the localized constructed analogs (LOCA) downscaling for the General Fluid Dynamics Lab Coupled Model (GFDL-CM3) and National Center for Atmospheric Research Community Earth System Model (CCSM4) global climate models (GCMS) (Pierce, Cayan, and Thrasher 2014; Gent et al. 2011; Griffies et al. 2011). LOCA downscaled data were retrieved from the downscaled CMIP3 and CMIP5 Climate and Hydrology Projections archive (http://gdo-dcp.ucllnl.org/downscaled_cmip_projections/).

Performance of GCMs in the North American Great Lake region varies due to the mostly unaccounted for regional climate impact of the Great Lakes (Notaro, Bennington, and Vavrus 2015; Rood and Briley 2018). GCM model selection was guided by performance and sensitivity (magnitude of change under a given emissions scenario relative to other GCMs) of the models in the Great Lakes region (Byun and Hamlet 2018). Both GFDL-CM3 and CCSM4 were found to be well performing models for the Great Lakes region, and CCSM4 and GFDL-CM3 represent a less and more sensitive models, respectively. Inspection of LOCA daily simulations of GFDL-CM3 and CCSM4 for the 1980-2009 period show excellent alignment with observed conditions during the same period (Figure 3). There appreas to be fewer dry days in the LOCA dataset relative to the observed dataset. Monthly and seasonal precipitation totals showed good agreement between the LOCA and observed datasets and no bias-correction was performed.

Figure 3: Comparison of LOCA daily projections and observed values for the climate normal period (1980-2009) from the GFDL-CM3 (A-C) and CCSM4 (D-F) global climate models. Points show median value, bold and thin bar show 1 and 2 standard devations respectively, and filled polygons show distribution of values.

Stochastic Weather Generator

The objective of projecting wetland conditions under future climate scenarios is to understand the expected response, and the range of possible responses. The 30-year periods of daily projections provide a representation of the normal climate variation under each scenario (Baddour and Kontongomde 2007). Quantifying the range of weather patterns under future climate scenarios and the corresponding wetland responses is best performed with many more years of simulation using stochastic processes (Daniel S. Wilks 2012). Stochastic weather generators (SWGs) are a tool for generating synthetic time series of weather that simulate conditions under an observed or projected climatology (D. S. Wilks and Wilby 1999). SWGs have also been used as a downscaling technique for global or regional climate models (Verdin et al. 2015, 2019). We followed the generalized linear model (GLM) based SWG described in Verdin et al. (2015), which is a type of the more general Richardson SWG (Richardson 1981).

Our SWG consists of four GLM models to simulate precipitation occurrence, precipitation amount, minimum daily temperature, and maximum daily temperature. Minimum and maximum daily temperatures are fit using separate models with the same model form. Each daily value is predicted using the previous day’s minimum and maximum temperatures, a pair of Fourier harmonics, and seasonal mean minimum and maximum temperatures. To simulate the daily variation inherent within a climate season additional error is added to the predicted daily minimum and maximum temperatures. The additional error is drawn from a multivariate skew-normal distribution conditioned on precipitation occurrence. This approach captures the correlation between daily minimum and maximum temperatures and between daily temperatres and precipitation occurrence. Precipitation occurrence is modeled using the logistic family with a probit link function. Occurrence is considered as a first order Markov process dependent upon the previous day’s precipitation occurrence (D. S. Wilks and Wilby 1999). Occurrence models also include a pair of Fourier harmonics to capture any seasonality in precipitation patterns, and the seasonal totals of precipitation for the observed period as covariates. Precipitation amounts are drawn from a gamma distribution defined by shape and scale parameters. The shape parameter is extracted from a gamma-family GLM with a log link fit to the observed precipitation amounts. The predictors in the model are a pair of Fourier harmonics and the seasonal precipitation totals. The scale of gamma distribution describing the distribution of storm sizes is a function of the GLM gamma-family shape and the predictors from GLM fit to precipitation amounts. This approach allows the storm distribution to vary temporally and capture seasonality in not only total precipitation but also the size of the storms. When seasonal means and totals are included as predictors within the GLM predictors it is said that the SWG is conditioned on the observed climate. When applying a conditioned SWG to a new climate projection it is assumed that the distribution of the precipitation amounts will remain constant under the new seasonally-defined climate scenario (Daniel S. Wilks 2012). We removed this assumption from our analysis by developing separate SWGs for each model-scenario combination using the LOCA daily projections for GLM fitting.

The work in Verdin et al. (2015) shows that these SWGs can capture and simulate the spatial correlation in temperature and precipitation. We did not include spatial correlation in our SWGs because of the relatively small size of the study area and the fact that each wetland is considered an independent system. For each model-scenario combination the SWG was fit using LOCA daily projections from GHCN-D sites USC00200718, USC00206220, USC00204104, USC00476398, USC00207812, USC00201088, USC00204328, USC00475352, USC00206215, USC00476398, USC00476518, USC00476939, USC00477140, USC00208680 to calculate annual regional seasonal precipitation and temperature for conditioning data. The projected daily data from the coordinates of the Bergland Dam meteorological station were used to fit the SWG GLMs. Each SWG was used to simulate 10,000 individual synthetic weather series with seasonal conditioning drawn randomly from the 30 years of projection data. Conditioning the models on individual years of observed data is intended to increase interannual (or inter-simulation) variability, which can otherwise be limited in SWGs (D. S. Wilks and Wilby 1999).

Data Analysis

Tests for statistical significance are unperformed and unreported with a single exception. This research is entirely dependent on simulation data drawing conclusions from modeled wetland water levels driven by synthetic weather series generated from parametric descriptions of LOCA downscaled of GCM projections. Reporting statistical significance and p-values would provide a false sense of certainty and dichomtomy to the future climate and wetland conditions. In our results we report and account for the uncertainty inherent in this type of simulation work. We draw conclusions by contrasting the probability of observing certain wetland condutions under different vegetation types and climate scenarios.

Total, EAB, and Climate Impacts

The distribution of modeled wetland water levels was evaluated for each future vegetation and climate scenario combination. The 10,000 simulations for each combination were summarized by simulated day of year to compute the median and the 67% highest density intervals (McElreath 2020) for each combination. This approach was used to quantify future conditions and baseline scenarios for comparison and benchmarking. Modeled baselines provide consistency between the baseline, or control, period and alternative vegetation cover and climate scenarios, which is critical to drawing meaningful conclusions. Apart from the advanatage of consistency, modeled baselines can also provide a flexible tool to answer more questions about the drivers of the the observed changes. Different baselines can be computed from the simulations by combining the six vegetation-climate combinations. To evaluate total impact of EAB and climate change, future non-black ash conditions were compared to black ash under the current climate. Alternative baseline comparisons include comparing each vegetation cover to itself under alternative climate scenarios, and comparing vegetation covers to each other within a climate scenario. These two alternative baselines allowed us to attribute how much of the observed total impact was attributable to each driver.

Critical Ecohydrological Thresholds

Observed water levels within these wetlands show high interannual variation under field control and treatment conditions. Rather than consider this low-signal, high-variance variable, we defined critical ecohydrological thresholds (CEHTs) as a measure of impact. CEHTs were set to capture wetland connectivity to the downstream hydrologic network via streamflow and subsurface flow, inundation when wetland water levels were near the soil surface with surface water likely in microtopography hollows, and drawdown when wetland water levels dropped far below the surface of the wetland (??). Wetland water levels were compared to these thresholds and the number of days that a wetland was above or a below a given threshold could be used to calculate the probability of occurrence of that CEHT. We chose to calculate these probabilities at the monthly scale, though can than be computed for time scales from daily to annually.

Results

SWG Performance

Synthetic weather generator performance can be evaluated by comparing generated weather sequences with the underlying observed/simulated climate data. The summary statistics of the sequences should align well with underlying data and the distribution of observed values should closely match the shape of the underlying data (Gregory, Wigley, and Jones 1993). Table 3 shows seasonal and climate model scenario summaries of LOCA and SWG variables. Values are reported as the mean and a 95% confidence interval (1.96 times standard deviations) around the mean. The percent of outliers were calculated using the intervals defined by the LOCA dataset. When comparing all of the simulated SWG data (equivalent to 10,000 years) the mean values and interquartile ranges of daily minimum and maximum temperature and precipitation align well with the LOCA-generated daily values (Table 3). Fewer total monthly precipitation outliers were observed for the SWG synthetic series than the LOCA climatology. Minimum and maximum temperatures showed an increase in outliers in the synthetic weather series for JJA and SON seasons, and a decrease in outliers for the DJF and MAM seasons. The difference in outliers was small and may be partially explained by the scale of the two datasets. The LOCA data represents 30 years of data compared to the 10,000 years of synthetic weather series. To compare the SWG and LOCA values as climatologies, we randomly sampled 30 sets of 30 years from the observed 10,000 synthetic series. The distribution of these samples were compared to the distribution of LOCA values in Figure 4. Figure 4A shows good agreement for all seasons for daily minimum and maximum temperature. The distribution of simulated monthly precipitation totals are shown in Figure 4B. Relative to minimum and maximum temperature, monthly precipitation totals show greater variability among the 30 sets of 30-year synthetic weather series. All three variables show that the synthetic weather climatologies cluster around the LOCA value as a mean.

Figure 4: Density plots of 30 years of LOCA generated data (bold), and 30 random samples of 30-year periods from the full 10,000 years of SWG synthetic weather series (light). Bold lines represent LOCA climatology, and fine lines show individual 30-year samples of SWG data.

Wetland Water Level Model Performance

Wetland water level model performance was evaluated using the retained independent testing datasets. We calculated R² as a metric for the relationship between observed and modeled values. R² cannot be used to evaluate the accuracy of the models because consistent over or under predictions can still result in high R² values (Krause, Boyle, and Bäse 2005). We use the median prediction error to measure model bias and the root median squared error (RMedSE) to assess overall predictive accuracy. RMedSE was used in place of RMSE becuase we expect some outlier errors where regionally-informed rainfall records do not match observed wetland rainfall. In addition to the overall accuracy of the models we want the models to accurately capture the probability of occurrence of the CEHTs. We tested the predicted probabilities of inundation, connectivity, and drawdown to the observed probabilities of the same conditions.

One wetland, 119, had all of the withheld training years with R² below 0.6, which was chosen as a threshold for unsatisfactory model results based on work by Moriasi et al. (2015) (Figure 5). Models at this wetland did not show noticeably higher rates of bias or error relative to the other wetlands. Site 119 is a closed basin with no surface outlet and was shown to have different source water characteristics than other sites (Van Grinsven et al. 2017). Table 4 summarizes the wetland model metrics by site status, presenting the minimum, maximum, mean, and median values observed across all site-year combinations. Mean and median model performance was the similar between site statuses for all model metrics (4). After excluding wetland 119 R² values ranged from 0.44 0.95.

Figure 5: Wetland model performance metrics. Metrics were calculated for each site-year combination and are presented within sites. RMedSE and rRMedSE are root median squared error and root median squared error relative to wetland water level range within that year.

When compared to observed wetland water levels, modeled wetland water levels showed no significant difference in probability of occurrence of CEHTs (Table ??). Reported values and significant tests are from a linear mixed model with dependent variable probability of occurrence, population effects of observed/predicted and wetland status, and a group effect for wetland (Lenth 2021). Summary and statistical tests were performed using estimated marginal means (Bates et al. 2015). The range of site-level probability of each level of interest were similar between modeled and observed data (Figure ??). Though non-significant these results show a slight systematic bias towards drier conditions in the modeled wetland water levels. All remaining comparisons of current conditions to combinations of future vegetation conditions and climate scenarios will use modeled current conditions. This comparison ensures that the described non-significant model bias does not impact results or conclusions. It also sets equal sample sizes between the control and future conditions.

Try vertical text for variable label then can horizontally expand LOCA and SWG values

Figure 6: Tests comparing predicted probability of occurrence for critical ecohydrological thresholds (CEHT). Connectivity, Inundation, and Drawdown respectively represent wetland water levels equal to or greater than modeled maximum water level, water levels within 10 cm of the wetland surface, and water levels more than 50 cm below the wetland surface. Reported values show the seasonal probability of daily water level exceeding the respective CEHT. Control and treated conditions refer to observed control and treated conditions or predicted modeled black ash and modeled non-forested conditions. Black point and error bar represents estimated marginal mean and its 95% confidence interval. Individual points represents data observed or simulated for 1 year at 1 site. All data presented derived from only the withheld test period for each wetland.

Wetland Water Levels

Wetland water levels under future climate conditions were highly variable under both scenarios and both vegetation types (Figure 7). Figure 7 shows the probability that wetland water levels will be above or below the CEHTs defined above. The reported probability is the proportion of observations where daily water level surpasses each threshold out of all modeled daily water levels in a month (~240,000). The median probability and a range representing 67% of all observations (highest density continuous intervals, HDCI) are shown as point and line estimates for the less and more sensitive future climate scenarios. Modeled black ash conditions are also reported as the median and 67% HDCI, represented as a crossbar and shaded area. Figures 10 & 9 present other comparisons using the same symbology structure.

For non-forested conditions both future climate scenarios showed a decrease in probability of connectivity relative to current wetland conditions (Figure 7A). The less sensitive climate scenario consistently showed a lower median probability of connectivity and was more likely to have simulated conditions outside of the HDCI of current conditions. The more-sensitive climate scenario showed greatest agreement with current conditions in July and August. Future forested conditions were more likely to show connectivity in July and August under the more senstive climate scenario, and only slightly less likely than current conditions under the less sensitive climate scenario (Figure 7B). Connectivity to the larger hydrologic network was much less likely to occur under the more sensitive climate scenario and both vegetation conditions in May (Figure 7A & B). Under both vegetation conditions the less-sensitive climate scenario resulted in a lower probability of connectivity relative to the more-sensitive climate scenario. While both vegetation conditions showed significant overlap between the less-sensitive climate conditions and control conditions, the non-forested conditions tended to skew towards current drier years (lower probability of connectivity) and future forested conditions skewed towards current wetter years (high probability of connectivity).

Inundation results in Figure 7C & D closely resemble the patterns observed in the connectivity results. Probability of inundation and connectivity are not exclusive. This is inline with fit model parameters that showed maximum sustained water levels were above the surface for almost every wetland and treatment combination (TABLE S1). Generally, inundated conditions would be less prevalent under future non-forested conditions than the future-forested conditions. For non-forested conditions there was little overlap between baseline conditions and the less sensitive climate scenario results (Figure 7D). This was contrary to non-forested conditions where the less-sensitive climate scenario results showed more overlap with baseline conditions, and the more senstive scenario varied by month relative to baseline conditions. Future-forested conditions showed increased probability for inundation in July, August, and September and both vegetation conditions showed lower probability for inundation in June (Figure 7D).

The probability that water levels will drop to less than 50 cm below the wetland soil surface agree with the connectivity and inundation results (Figure 7E & F). Under both climate scenarios the future forest conditions showed a lower probability of drawdown than under baseline condtiions, with almost no overlap with current conditions in July, August, and September (Figure 7F). Drawdown events were very unlikely with near zero probability of occurring under future forested conditions under both climate scenarios. Probability of drawdown in non-forested conditions increased in magnitude and variability dramatically from July through August under both future climate scenarios. Differences from current conditions were greatest in August under the less senstive climate scenario (Figure 7E).

Figure 7: Probability of observing water levels above (connectivity, inundation) or below (drawdown) ecohydrologically significant thresholds. Points represent median the median observed probability within a month and ranges represent 67% of the simulations as highest-density continuous interval (HDCI). Gray shaded area represents model simulations for control black (ash forest) conditions under the current climate and the black bar represents the control monthly median probability.

Discussion

Talk about more senstive connectivity in May and inundation in June as indicators of earlier melt Can tie it back to the sourcewater work for when snowmelt signal was present

SWG

Our SWG performed well at simulating weather series that conformed to LOCA-downscaled future climate scenarios (Table 3, Figure 4). The synthetic weather series were well aligned with the LOCA series, with precipitation showing more inter-simulation variability than minimum or maximum temperature. This result is understandable when we consider that in the study region precipitation is only weakly seasonal, while minimum and maximum temperature have strong seasonal signals. There may be a slight trend in under estimating the range of observed future climate conditions (Table 3). This is likely a result of reversion to the seasonal mean due to a lack of available external weather drivers. Under future climate scenarios extreme events are expected to be more commonplace (CITATION), and these results suggest event size could be better parameterized by increasing the probability of extreme events through selection of distribution choices and modeling framework (Verdin et al. 2019). This shortcoming of the SWGs should not affect our findings because our objective is not to model response to extreme events, but overall probability of CEHT water levels. These water levels will be affected by extreme precipitation events but the systems do not have infinite storage and more extreme events can be expected to pulse through the system more quickly (CITATION). The quicker the pulse passes through the system the smaller the impact on overall probability of wetland water levels exceeding CEHTs

We chose to use an SWG approach that relied on only a single estimate of climate within the study area. Alternatively, regional variation in climate and weather patterns could be incorporated by using additional meteorological stations in the development of the SWG. A spatially-explicit SWG approach takes into account the covariance of regional meteorological stations and has been shown to faithfully capture larger regional trends in climate and weather (Verdin et al. 2015). Within our study region there are no clear weather or climate gradients and so we assumed that each study wetland was within a single climatic unit. This results in a much simpler SWG-construction process, but the resulting SWG cannot be used outside of the study region without reconstruction from a different meteorological station(s). To simulate black ash wetland response on a regional basis would require addressing known temperature and precipitation gradients through the range of black ash.

In addition to using only a single spatial representation of future climate, we chose to use individual years of simulation. The direct impact of this is to missing long-term droughts and wet spells. These events are expected to increase in frequency under future climate scenarios (CITATION). We deliberately chose to use single-year simulations because of the length of available training data. During our study the region was in the process of transitioning from a relatively dry period to a relatively wetter period (@ref{fig:water-balacne}). Our training data are primarily drawn from the drier portion of the study period. Without longer-periods of training data we do not feel that we can capture both the intra- and inter-annual dynamics of wetland hydrology. If we assumed all 10000 years of synthetic weather series were contiguous weather records, the longest period of drought years (annual precipitation less than annual PET) was 7 and 9 years for the less and more sensitive climate scenarios, respectively. By using single-year simulations we are implicitly assuming that fall and winter (SON and DJF) precipitation is high enough to recharge the local hydrology driving these wetlands, restarting the cycle at or near maximum wetland water levels. When we consider only fall and winter precipitation we found 0%, 8.7%, and 10.1% of synthetic weather series had totals below the tenth percentile of observed data, less-sensitive, and more-sensitive climate scenarios, respectively. Unsurprisingly close to 10% of the synthetic weather winters were drier than the tenth percentile of the LOCA downscaled values. Synthetic (and LOCA) winters were wetter than observed winters, supporting the assumption that dormant-season precipitation would be enough to recharge wetland water levels.

Figure 8: Regional cumulative water availability for the period from November 1, 2005 through October 21, 2020. Water availability was calculated as the difference between rainfall plus melt and potential evapotranspiration.

Ecosystem Specific Yield

The approach used to calculate E_Sy for use in relating the magnitude of water level change to the drivers of that change was an alternative to those presented in previous research. We developed this alternative based on unsuccessful attempts to implement E_Sy calculates through the rainfall/rise method described in McLaughlin, et al (-McLaughlin and Cohen (2014)). Directly proving the accuracy of this approach is outside the scope of this study, and would be best performed through theoretical and designed experimental work. Our research has however provided empirical evidence for the application of our approach. Section CASE_STUDY in Chapter TRANSDUCER_ERRORS demonstrates that ET estimates derived using this approach perform comparably to other implementations for calculating E_Sy (see White method details above). Those results showed that calculated ET was correlated with daily PET estimates from nearby meteorological stations. In this work we can look at how well a simple regression between daily water level change and adjusted or unadjusted hydrologic drivers perform. Figure ?? shows the results of a quantile regression (not used elsewhere in this study) explaining daily change in wetland water level as a function of PET P and M. We found that R² of the model improved from 0.69 to 0.78. While 0.69 is considered an acceptable R² as set out in Moriasi (-Moriasi et al. (2015)), we can see that water level drawdown is extremely underestimated by the unadjusted model (Figure ??A). Our results show that periods with the highest drawdown are not correctly modeled using the raw inputs. These periods of high drawdown are captured when the E_Sy adjusted hydrologic inputs are used. Periods of high drawdown generally occurred mid-season when water levels were farther below the surface. Our approach agrees with others (CITATIONS) finding that the impact of E_Sy increases as water levels drawdown and that flooded conditions approach an E_Sy of 1 (CITATION). This relationship results in the largest adjustment to inputs during these mid-season periods, explaining the differences seen between panels A and B in Figure ??. These results together with those in Section CASE_STUDY in Chapter TRANSDUCER_ERRORS provide empirical evidence that this alternative formulation captures the variation of E_Sy with wetlands and can be used in scenarios where the rainfall-runoff ratio approach is masked by outside influences such as low-frequency hydrologic signals (Zhu et al. (2011)) and surface water connectivity (Watras et al. (2017)).

What if the response of black ash transpiration to lowering water levels (rapid increase in T) is actually the signal I am capturing, not E_Sy? This may be partly hapeening, but the response between rainfall and water level response is well explained by the current E_Sy response, suggesting that the modeled phenomenon is primarily E_Sy.

Wetland Water Level Models

Wetland model performance was adequate for the objective of this study, which was to observe changes in the hydrographs of black ash wetlands and compare the probability of critical ecohydrological water levels occurring. Overall the retained models showed good correlation between observed and modeled water levels (Figure 5A, Table @ref(tab:model_metrics)). We found that 80% and 100% of site-year combinations outside of sites 119 and 156 had R^2 values above 0.6 for the control and treated condition models, respectively. Model performance varied by site with a median rRMedSE of 12.75% (Control: 13.03%, Treated: 12.58%) of annual water level range and median model error was -2.19 cm (Control: -1.93 cm, Treated: -2.28 cm). RMedSE and relative RMedSE had similar patterns among sites (Figure 5). This suggests that from both absolute and relative standpoint, the wetland models as developed capture dynamic wetland water levels more accurately than static wetland water levels. Daily percent bias estimates (PBIAS) would be helpful to evaluate error relative to observed daily water level. However, wetland water levels cross 0 so that the relative error approaches infinity as water levels rise or fall towards the surface.

The wetland models performed well predicting the probability of occurrence of important hydroecological water levels. Table 5 and Figure ?? show how probability of occurrence of connectivity, inundation, and drawdown were not found to be different between the observed and modeled data in the test periods. They also show that there is no non-significant systematic model bias (positive or negative). Although no systematic model bias was identified, use a modeled control baseline is still an important safeguard against drawing conclusions based on potential from model artefacts. Comparing the systematically-biased results could exaggerate or mask real expectation of wetter or drier wetland conditions.

The next step for these models is to fit them using a hierarchical Bayesian approach. This should improve model fit by incorporating information from all sites into a single model for each vegetation condition while still allowing model parameters to vary between sites. It will also increase the available data for model fitting as unbalanced vegetation-conditions would be accounted for in the structure of the model. If model performance does not significantly improve with that approach more, or potentially all years of data can be used to parameterize the wetland models. Using most or all of the data for model parameterization precludes the options of evaluating model performance on truly independent data. Kozak and Kozak (-Kozak and Kozak (2003)) provide an alternative view that replaces cross-validation with lack-of-fit metrics, which evaluate model predictions for assumptions of model assumptions. Their main conclusion was that cross-validation increase the variance of model estimates potentially resulting in poor model choice.

Future Hydrology

The importance of these simulation results is to help relate potential future simulations to current conditions in dry, wet, or ‘normal’ years. This allows researchers and managers to answer questions similar to: What will these sites look like in a ‘normal’ year in the future? How will these sites respond to wet/dry conditions in the future? For example, we could expect that during a wet year under the more sensitive future climate scenario a non-forested wetland on these sites would have a similar probability of surface inundation as a black ash wetland does today (Figure 7C). In general, black ash wetlands that become non-forested can be expected to be have wetter conditions under less- sensitive future climate scenarios, and under more-sensitive future climate scenarios are more likely to experience conditions similar to today. Conclusions could be reached for each combination of climate and vegetaiton scenario and this information can and should be used to influence management approaches and timing in responding to EAB in black ash wetlands.

The combined impact of EAB and climate can be seen by comparing the future climate conditions (blue and orange) to the control conditions under historical climate (shaded gray) in Figure 7. We can draw additional information about the individual effects of future climate and vegetative cover on wetlands by altering how we calculate the baseline. In Figure 7 the baseline is to simulated black ash forests under current climate conditions. To isolate the effects of each climate scenario of wetland hydrology we can compare the response of a vegetative under future climate conditions against its response under current climate conditions. Conversely we can evaluate the impact of vegetative cover by comparing each vegetative cover to simulated black ash forests under current climate conditions. Figures 9 and 10 show the probability of occurrence of each critical ecohydrological water level under conditions where the vegetative cover and climate, respectively, are held constant. Table 6 presents the change in probability of occurrence for each critical water level and alternative vegetative condition that can be attributed to EAB or a less or more sensitive future climate scenario.

We found that EAB impacts lead to slightly wetter conditions under non-forested conditions and the present climate. This agrees with observed water level response analysis in Van Grinsven, et al (-Van Grinsven et al. (2017)), where absolute water level response to simulated post-EAB conditions resulted in wetter conditions and significantly lower growing season drawdown rates. A modeled alternate forest composition (Future Forested) shows that EAB impact alone would lead to dramatically wetter conditions relative to black ash forests. This result is expected as we modeled our future forested conditions based on previous work showing that co-dominant hardwoods had significantly lower seasonal transpiration estimates than black ash (Shannon et al. (2017)). Under both alternative vegetative conditions we found that the isolated impact of future climate scenarios would lead to much drier sites. The probability of surface water connectivity or inundation occurring are dramatically decreased under both climate scenarios. The probability of wetland water levels dropping below 50 cm were increased under future climate scenarios for the non-forested conditions, but showed little change for non-forested conditions.

Figure 9: Probability of occurrence for critical ecohydrological wetland water levels under combinations of vegetative cover and climate conditions. Each column compares simulated wetland response within a fixed vegetation condition across multiple climate scenarios. Probability of occurence is the proportion of days in each month simulated water levels (10000 simulations) reached each threshold water level. Connectivity: connected to surface drainage; Inundation: within 10 cm of soil surface; Drawdown: more than 50 cm below soil surface. Values are reported as the median and the bounds of a 67% highest density continuous interval (an interval, potentially asymmetric, that contains 67% of simulated values).

Figure 10: Probability of occurrence for critical ecohydrological wetland water levels under combinations of vegetative cover and climate conditions. Each column compares simulated wetland response within a fixed climate scenario across multiple potential vegetative covers. Differences within a column represent the impact of EAB and potential management decisions. Probability of occurence is the proportion of days in each month simulated water levels (10000 simulations) reached each threshold water level. Connectivity: connected to surface drainage; Inundation: within 10 cm of soil surface; Drawdown: more than 50 cm below soil surface. Values are reported as the median and the bounds of a 67% highest density continuous interval (an interval, potentially asymmetric, that contains 67% of simulated values).

Our hypothesis that the impacts of EAB and climate change would counteract each other is partially supported by our hypothesis. The dynamics of wetland water levels are complex (CITATION) climate change and invasive species can each have major consequences (CITATIONS). We can summarize our findings by saying we expect wetlands that transition from black ash forests to non-forested sites to become drier, primarily due to the effects of increased summer evaporative demand (Figures 7 & 11, Table 6). Our simulations saw less summer precipitation under the less sensitive scenario, while under the more sensitive scenario we obseved an increase in large storms (Figure 11. If management or natural regeneration leads to establishment of future forests with lower site transpiration rates, these wetlands can be expected to have the occurrence of high water conditions remain stable or increase in frequency, while dry conditions (water levels more than 50 cm below the surface) become more rare (Figures 7 & 11, Table 6). Our results show that the impact of climate change in this region will lead to consistently drier wetland conditions. Changes to vegetative cover can amplify or counteract hydrologic alterations attributable to climate change.

Figure 11: Distribution of 30-day total precipitation and potential evapotranspiration presented by climate season and climate scenario. All data present is fro the Bergland Dam, MI, USA GHCND station. Observed climate is from 1980-2009, future climate is from LOCA-downscaled climate scenarios for the period from 2070-2099. Potential evapostranspiration was calculated from observed and GCM outputs using the Hargreaves-Samani equation ((Hargreaves and Allen 2003)).

Drivers of Future Wetland Hydrology

LOCA-downscaled climate simulations for the study region showed

It is expected that vegetation conditions with lower transpiration rates have lower rates of water level drawdown and therefore lower increased water levels. In these wetlands this effect is compounded by the dynamics of E_Sy. E_Sy increases in magnitude as water levels decline (Figure @(fig:esy-function)), which may be a result of reduced soil pore space and wetland cross-sectional area so that a smaller volume of change results in a large water level difference (CITATIONS). The effect is that PET results in smaller changes to wetland water levels under wet conditions than dry conditions. In this way the impact of EAB may lead to a feedback loop in which water levels remain elevated:

1. AET begins to draw down water levels, E_Sy is low,
   
2. Water levels decline slowly because of reduced AET and low E_Sy,
   
3. E_Sy is lower than under black ash as PET increases to mid-season peak,
   
4. High PET impact on water levels are reduced by lower AET and sustained low E_Sy.

Under the current black ash conditions E_Sy increases more rapidly due to higher black ash AET. The increase in E_Sy accelerates the impact of higher AET, creating faster and larger declines in wetland water level throughout the growing season (Figure ??). This feedback loop may partially explain the persistence of hydrologic impact following EAB disturbance observed in Michigan and Minnesota (Diamond et al. 2018; Kolka et al. 2018). In effect, the impact of EAB may shock these systems into an alternative stable state of elevation water levels (Scheffer and Carpenter 2003). Our results indicate that future climate scenarios will likely have a large enough impact to again shock these systems out of a stable state. We cannot say from these simulations whether the wetlands will reach a new stable state under future climate conditions or what that state would be.

Make sure talking about PET and P of future climates vs current climate in reference to future-climate figure

We observed wetter conditions observed under our Future Forest simulations relative to both Black Ash and Non-Forested simulations under all climate scenarios (Figures ??, 10). Two factors are likely contributing to this result. The first is that our future forest composition is known to have lower transpiration rates than the existing black ash canopy (Shannon et al. 2018). As a model choice we opted to assume that the next most dominant canopy species would be the likely replacement canopy species. The future forest we modeled is a narrow range of potential future forest compositions. natural regeneration or planting efforts may lead to future forested species composition that more closely matches the evaporative potential of the current black ash canopy. Secondly, the future forested conditions may have a lower total AET than the non-forested conditions. In black ash wetlands in Minnesota simulated post-EAB conditions (girdled and standing ash) were found to have higher water levels than site were the ash stems were harvested and removed (Diamond et al. 2018). The authors attributed the result to limited AET leading to reduced solar energy and wind-driven boundary layer mixing due to the still standing stems. Our conditions have a notable difference from that study in that we are assuming there are living trees on site. However a closed canopy would have the same effect of reducing open water evaporation and understory transpiration. With suffificiently reduced canopy transpiration the effect of certain forest compositions could result in reduced overall AET. This was observed in reverse in northern Minnesota when total AET increased after low-transpiration black spruce was removed from a peatland (CITATION).

The contrast between the Non-Forested and Future-Forested simulated hydrology suggests an important management tool. Iverson et al (2016) and the work from Bolton (2018) and Looney (2015) laid out a framework and results for evaluating potential replacement species considering site conditions of black ash wetlands. The results presented here show the opportunity for management decisions that take into account the impact of future vegetation on site hydrologic conditions. The general trend of drier future conditions can be to some degree counteracted by management for cover with lower evapotranspiration rates. This tactic could be used to retain water on the landscape, creating refugia of standing water or cool moist soils for flora and fauna. Although much of the region is expected to have drier summer our localized study area is expected to have wetter summers. Our results can still be applied to the broader region. The impact of climate change in other areas of the region would become even more pronounced. Vegetation conditions would still amplify or counteract climate impacts, and may become even more important to providing wet or moist refugia on the landscape.

Future Research

It should be stated that these models and simulations cannot capture the full range of potential changes from EAB, other potential invasive species, and climate change. The following is a non-exhaustive list of further known assumptions within this analysis that could warrant further study.

• Watras et al. (2014) found that there are decadal scale oscillations of water levels within the Great Lakes and inland lakes in the Great Lakes Region. These can be expected to impact wetland water levels similarly.
  
• While our results showed a smaller and earlier snowmelt pulse (TODO: output figure) we did not attempt to model the more complex dynamics of rain-on-snow events and a compressed melting period reducing recharge to the intermediate groundwater sources shown to feed these wetlands (snowmetl citation, GROUNDWATER CITATION Van Grinsven et al. 2017).
  
• Our single-year simulations do not account for the potential effects of multi-year drought or for periods of significantly reduced fall/winter precipitation on wetland water level rebound.
  
• The response of each species and plant community will show individual climate change responses that are not captured in these models and may be non-linear (Short 2016).

Beyond the potential for unforeseen vegetation responses to future climates is the uncertainty of future climate conditions. The authors of the fourth national climate assessment have highlighted that early climate prediction have under-predicted contemporary shifts in response to climate change (CITATION). Unfortunately this means that our models built on those simulations may be underestimating the magnitude of future changes.

Wetlands provide a range of ecological services from water quality and stormflow retention to carbon storage, all of which stand to be impacted by climate change (Moomaw et al. 2018). We did not attempt to quantify how changes in wetland hydrology will influence these other services. By highlighting the magnitude and drivers of hydrologic change in these systems we hope to spur future research.

Conclusions

“3) site conditions following EAB infestation, and 4) site conditions in a future climate.”

Our research has shown that changes in evaporative demand and precipitation regimes will likely result in drier conditions on what are now black ash wetlands. In addition the functional loss of black ash due to EAB presents challenges and opportunities for the future of these wetlands. A large extent of forested wetlands may require intervention to retain desired benefits or features they currently provide. These interventions can be used to drive sites towards wetter conditions, retaining water on the landscape that may otherwise be lost under a changing climate.

To implement long-term management objectives in wetland and forest management requires a deep understanding of the systems. We have shown that as water levels decline, the impact of each additional driver increases due to changes in E_Sy. Changes in AET can interact with E_Sy to create feedback loops underscoring that knowledge of species adaptation to wet conditions and capacity to respond quickly to drying conditions is critical for projected future site conditions. The transition from black ash to alternative vegetative cover can counteract or amplify the impacts of future conditions with higher evaporative demand. And the relationship between evaporative demand and E_Sy can amplify the effect of wetland water level drawdown.

Acknowledgements

I would like to thank Andrew Verdin for discussion of synthetic weather generators and code examples. Additionally my wife Danielle Shannon for both the support she provided during writing, and for share her knowledge of evaluting and communicating climate change impacts on forests.

References

Arguez, Anthony, Imke Durre, Scott Applequist, Russell S. Vose, Michael F. Squires, Xungang Yin, Richard R. Heim, and Timothy W. Owen. 2012. “NOAA’s 1981 U.S. Climate Normals: An Overview.” Bulletin of the American Meteorological Society 93 (11): 1687–97. https://doi.org/10.1175/BAMS-D-11-00197.1.

Baddour, Omar, and Hama Kontongomde, eds. 2007. “THE ROLE OF CLIMATOLOGICAL NORMALS IN A CHANGING CLIMATE.” WCDMP-No. 61 , WMO-TD No. 1377. Geneva: World Meteorological Organization.

Bates, Douglas, Martin Mächler, Ben Bolker, and Steve Walker. 2015. “Fitting Linear Mixed-Effects Models Using lme4.” Journal of Statistical Software 67 (1): 1–48. https://doi.org/10.18637/jss.v067.i01.

Bojanowski, Jedrzej S. 2016. “Sirad: Functions for Calculating Daily Solar Radiation and Evapotranspiration.” https://CRAN.R-project.org/package=sirad.

Bolton, Nicholas, Joseph Shannon, Joshua Davis, Matthew Grinsven, Nam Noh, Shon Schooler, Randall Kolka, Thomas Pypker, and Joseph Wagenbrenner. 2018. “Methods to Improve Survival and Growth of Planted Alternative Species Seedlings in Black Ash Ecosystems Threatened by Emerald Ash Borer.” Forests 9 (3): 146. https://doi.org/10.3390/f9030146.

Brinson, Mark M. 1993. “A Hydrogeomorphic Classification for Wetlands.” Wetlands Research Program Technical Report WRP-DE-4 WRP-DE-4 (August): 101. https://doi.org/10.2134/agronj2001.931131x.

Bristow, Keith L., and Gaylon S. Campbell. 1984. “On the Relationship Between Incoming Solar Radiation and Daily Maximum and Minimum Temperature.” Agricultural and Forest Meteorology 31 (2): 159–66. https://doi.org/10.1016/0168-1923(84)90017-0.

Burkett, Virginia, and Jon Kusler. 2000. “CLIMATE CHANGE: POTENTIAL IMPACTS AND INTERACTIONS IN WETLANDS OF THE UNITED STATES ¹.” JAWRA Journal of the American Water Resources Association 36 (2): 313–20. https://doi.org/10.1111/j.1752-1688.2000.tb04270.x.

Bürkner, Paul-Christian. 2017. “brms: An R Package for Bayesian Multilevel Models Using Stan.” Journal of Statistical Software 80 (1): 1–28. https://doi.org/10.18637/jss.v080.i01.

———. 2018. “Advanced Bayesian Multilevel Modeling with the R Package brms.” The R Journal 10 (1): 395–411. https://doi.org/10.32614/RJ-2018-017.

Byun, Kyuhyun, and Alan F. Hamlet. 2018. “Projected Changes in Future Climate over the Midwest and Great Lakes Region Using Downscaled Cmip5 Ensembles: PROJECTED CLIMATE CHANGES OVER THE MIDWEST AND GREAT LAKES REGION.” International Journal of Climatology 38 (April): e531–53. https://doi.org/10.1002/joc.5388.

Coron, L., G. Thirel, O. Delaigue, C. Perrin, and V. Andréassian. 2017. “The Suite of Lumped GR Hydrological Models in an R Package.” Environmental Modelling & Software 94 (August): 166–71. https://doi.org/10.1016/j.envsoft.2017.05.002.

Davis, Joshua C., Joseph P. Shannon, Nicholas W. Bolton, Randall K. Kolka, and Thomas G. Pypker. 2017. “Vegetation Responses to Simulated Emerald Ash Borer Infestation in Fraxinus Nigra Dominated Wetlands of Upper Michigan, USA.” Canadian Journal of Forest Research 47 (3): 319–30. https://doi.org/10.1139/cjfr-2016-0105.

Diamond, Jacob S., Daniel L. McLaughlin, Robert A. Slesak, Anthony W. D’Amato, and Brian J. Palik. 2018. “Forested Versus Herbaceous Wetlands: Can Management Mitigate Ecohydrologic Regime Shifts from Invasive Emerald Ash Borer?” Journal of Environmental Management 222 (September): 436–46. https://doi.org/10.1016/j.jenvman.2018.05.082.

Gent, Peter R., Gokhan Danabasoglu, Leo J. Donner, Marika M. Holland, Elizabeth C. Hunke, Steve R. Jayne, David M. Lawrence, et al. 2011. “The Community Climate System Model Version 4.” Journal of Climate 24 (19): 4973–91. https://doi.org/10.1175/2011JCLI4083.1.

Gregory, Jonathan M, TML Wigley, and PD Jones. 1993. “Application of Markov Models to Area-Average Daily Precipitation Series and Interannual Variability in Seasonal Totals.” Climate Dynamics 8 (6): 299–310.

Griffies, Stephen M., Michael Winton, Leo J. Donner, Larry W. Horowitz, Stephanie M. Downes, Riccardo Farneti, Anand Gnanadesikan, et al. 2011. “The GFDL Cm3 Coupled Climate Model: Characteristics of the Ocean and Sea Ice Simulations.” Journal of Climate 24 (13): 3520–44. https://doi.org/10.1175/2011JCLI3964.1.

Haack, Robert A, Eduard Jendek, Houping Liu, Kenneth R Marchant, Toby R Petrice, Therese M Poland, Hui Ye, and East Lansing. 2002. “The Emerald Ash Borer: A New Exotic Pest in North America,” 5.

Hargreaves, George H., and Richard G. Allen. 2003. “History and Evaluation of Hargreaves Evapotranspiration Equation.” Journal of Irrigation and Drainage Engineering 129 (1): 53–63. https://doi.org/10.1061/(ASCE)0733-9437(2003)129:1(53).

Hayhoe, Katharine, Jeff VanDorn, Thomas Croley, Nicole Schlegal, and Donald Wuebbles. 2010. “Regional Climate Change Projections for Chicago and the US Great Lakes.” Journal of Great Lakes Research 36 (January): 7–21. https://doi.org/10.1016/j.jglr.2010.03.012.

Herms, Daniel A., and Deborah G. McCullough. 2014. “Emerald Ash Borer Invasion of North America: History, Biology, Ecology, Impacts, and Management.” Annual Review of Entomology 59 (1): 13–30. https://doi.org/10.1146/annurev-ento-011613-162051.

Horel, John, Michael Splitt, L Dunn, J Pechmann, B White, C Ciliberti, S Lazarus, J Slemmer, D Zaff, and J Burks. 2002. “Mesowest: Cooperative Mesonets in the Western United States.” Bulletin of the American Meteorological Society 83 (2): 211–26.

“Hourly Precipitation Data (HPD) Network, Version 2.R2.” 2021. NOAA National Centers for Environmental Information.

Iverson, Louis, Kathleen S. Knight, Anantha Prasad, Daniel A. Herms, Stephen Matthews, Matthew Peters, Annemarie Smith, Diane M. Hartzler, Robert Long, and John Almendinger. 2016. “Potential Species Replacements for Black Ash (Fraxinus Nigra) at the Confluence of Two Threats: Emerald Ash Borer and a Changing Climate.” Ecosystems 19 (2): 248–70. https://doi.org/10.1007/s10021-015-9929-y.

Kolka, Randall, Anthony D’Amato, Joseph Wagenbrenner, Robert Slesak, Thomas Pypker, Melissa Youngquist, Alexis Grinde, and Brian Palik. 2018. “Review of Ecosystem Level Impacts of Emerald Ash Borer on Black Ash Wetlands: What Does the Future Hold?” Forests 9 (4): 179. https://doi.org/10.3390/f9040179.

Kozak, Antal, and Robert Kozak. 2003. “Does Cross Validation Provide Additional Information in the Evaluation of Regression Models?” 33: 12.

Krause, P., D. P. Boyle, and F. Bäse. 2005. “Comparison of Different Efficiency Criteria for Hydrological Model Assessment.” Advances in Geosciences 5 (December): 89–97. https://doi.org/10.5194/adgeo-5-89-2005.

Lenth, Russell V. 2021. Emmeans: Estimated Marginal Means, Aka Least-Squares Means. https://github.com/rvlenth/emmeans.

Loheide II, Steven P. 2008. “A Method for Estimating Subdaily Evapotranspiration of Shallow Groundwater Using Diurnal Water Table Fluctuations.” Ecohydrology 1 (1): 59–66. https://doi.org/10.1002/eco.7.

Looney, Christopher E., Anthony W. D’Amato, Shawn Fraver, Brian J. Palik, and Michael R. Reinikainen. 2016. “Examining the Influences of Tree-to-Tree Competition and Climate on Size-Growth Relationships in Hydric, Multi-Aged Fraxinus Nigra Stands.” Forest Ecology and Management 375 (September): 238–48. https://doi.org/10.1016/j.foreco.2016.05.050.

Looney, Christopher E., Anthony W. D’Amato, Brian J. Palik, and Robert A. Slesak. 2015. “Overstory Treatment and Planting Season Affect Survival of Replacement Tree Species in Emerald Ash Borer Threatened Fraxinus Nigra Forests in Minnesota, USA.” Canadian Journal of Forest Research 45 (12): 1728–38. https://doi.org/10.1139/cjfr-2015-0129.

Looney, Christopher E., Anthony W. D’Amato, Brian J. Palik, Robert A. Slesak, and Mitchell A. Slater. 2017. “The Response of Fraxinus Nigra Forest Ground-Layer Vegetation to Emulated Emerald Ash Borer Mortality and Management Strategies in Northern Minnesota, USA.” Forest Ecology and Management 389 (April): 352–63. https://doi.org/10.1016/j.foreco.2016.12.028.

McElreath, Richard. 2020. Statistical Rethinking: A Bayesian Course with Examples in r and Stan. 2nd ed. Boca Raton, Florida: Chapman; Hall/CRC.

McLaughlin, Daniel L., and Matthew J. Cohen. 2014. “Ecosystem Specific Yield for Estimating Evapotranspiration and Groundwater Exchange from Diel Surface Water Variation: ECOSYSTEM S _Y FOR ESTIMATING ET AND GROUNDWATER EXCHANGE.” Hydrological Processes 28 (3): 1495–1506. https://doi.org/10.1002/hyp.9672.

McLaughlin, Daniel L., Jacob S. Diamond, Carlos Quintero, James Heffernan, and Matthew J. Cohen. 2019. “Wetland Connectivity Thresholds and Flow Dynamics From Stage Measurements.” Water Resources Research 55 (7): 6018–32. https://doi.org/10.1029/2018WR024652.

Menne, Matthew J., Imke Durre, Bryant Korzeniewski, Shelley McNeal, Kristy Thomas, Xungang Yin, Steven Anthony, et al. 2012. “Global Historical Climatology Network - Daily (GHCN-Daily), Version 3.28-Upd-2021031419.” http://doi.org/10.7289/V5D21VHZ.

Menne, Matthew J., Imke Durre, Russell S. Vose, Byron E. Gleason, and Tamara G. Houston. 2012. “An Overview of the Global Historical Climatology Network-Daily Database.” Journal of Atmospheric and Oceanic Technology 29 (7): 897–910. https://doi.org/10.1175/JTECH-D-11-00103.1.

Moomaw, William R., G. L. Chmura, Gillian T. Davies, C. M. Finlayson, B. A. Middleton, Susan M. Natali, J. E. Perry, N. Roulet, and Ariana E. Sutton-Grier. 2018. “Wetlands In a Changing Climate: Science, Policy and Management.” Wetlands 38 (2): 183–205. https://doi.org/10.1007/s13157-018-1023-8.

Moriasi, D. N., M. W. Gitau, N. Pai, and P. Duggaupati. 2015. “Hydrologic and Water Quality Models: Performance Measures and Evaluation Criteria.” Transactions of the ASABE 58 (6): 1763–85. https://doi.org/10.13031/trans.58.10715.

Notaro, Michael, Val Bennington, and Steve Vavrus. 2015. “Dynamically Downscaled Projections of Lake-Effect Snow in the Great Lakes Basin.” Journal of Climate 28 (4): 1661–84. https://doi.org/10.1175/JCLI-D-14-00467.1.

Pierce, David W., Daniel R. Cayan, and Bridget L. Thrasher. 2014. “Statistical Downscaling Using Localized Constructed Analogs (LOCA).” Journal of Hydrometeorology 15 (6): 2558–85. https://doi.org/10.1175/JHM-D-14-0082.1.

R Core Team. 2019. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.

Richardson, C. W. 1981. “Stochastic Simulation of Daily Precipitation, Temperature, and Solar Radiation.” Water Resources Research 17 (1): 182–90. https://doi.org/10.1029/WR017i001p00182.

Rood, Richard, and Laura Briley. 2018. “Great Lakes Ensemble April 2017 - April 2018 Progress Report.” Great Lakes Integrated Sciences + Assessments. http://glisa.umich.edu/media/files/Ensemble_Progress_Report_April_2018.pdf.

Scheffer, Marten, and Stephen R. Carpenter. 2003. “Catastrophic Regime Shifts in Ecosystems: Linking Theory to Observation.” Trends in Ecology & Evolution 18 (12): 648–56. https://doi.org/10.1016/j.tree.2003.09.002.

Shannon, Joseph, Joshua Davis, Matthew Van Grinsven, Nicholas Bolton, Nam Jin Noh, Thomas Pypker, Randall Kolka, and Joseph Wagenbrenner. 2017. “Water Level Controls on Transpiration of Co-Dominant Species in Black Ash Wetlands.” Poster {{Presentation}}. Duluth, MN, USA.

Shannon, Joseph, Matthew Van Grinsven, Joshua Davis, Nicholas Bolton, Nam Noh, Thomas Pypker, and Randall Kolka. 2018. “Water Level Controls on Sap Flux of Canopy Species in Black Ash Wetlands.” Forests 9 (3): 147. https://doi.org/10.3390/f9030147.

Short, Frederick T. 2016. “Impacts of Climate Change on Submerged and Emergent Wetland Plants.” Aquatic Botany, 15.

Slesak, Robert A., Christian F. Lenhart, Kenneth N. Brooks, Anthony W. D’Amato, and Brian J. Palik. 2014. “Water Table Response to Harvesting and Simulated Emerald Ash Borer Mortality in Black Ash Wetlands in Minnesota, USA.” Canadian Journal of Forest Research 44 (8): 961–68. https://doi.org/10.1139/cjfr-2014-0111.

Swanston, Christopher W., Maria K. Janowiak, Leslie A. Brandt, Patricia R. Butler, Stephen D. Handler, P. Danielle Shannon, Abigail Derby Lewis, et al. 2016. “Forest Adaptation Resources: Climate Change Tools and Approaches for Land Managers. 2nd Ed.” NRS-GTR-87-2. Newtown Square, PA: U.S. Department of Agriculture, Forest Service, Northern Research Station. https://doi.org/10.2737/NRS-GTR-87-2.

Valéry, Audrey, Vazken Andréassian, and Charles Perrin. 2014. “‘As Simple as Possible but Not Simpler’: What Is Useful in a Temperature-Based Snow-Accounting Routine? Part 1 Comparison of Six Snow Accounting Routines on 380 Catchments.” Journal of Hydrology 517 (September): 1166–75. https://doi.org/10.1016/j.jhydrol.2014.04.059.

van Vuuren, Detlef P., Jae Edmonds, Mikiko Kainuma, Keywan Riahi, Allison Thomson, Kathy Hibbard, George C. Hurtt, et al. 2011. “The Representative Concentration Pathways: An Overview.” Climatic Change 109 (1-2): 5–31. https://doi.org/10.1007/s10584-011-0148-z.

Van Grinsven, Matthew J., Joseph P. Shannon, Joshua C. Davis, Nicholas W. Bolton, Joseph W. Wagenbrenner, Randall K. Kolka, and Thomas G. Pypker. 2017. “Source Water Contributions and Hydrologic Responses to Simulated Emerald Ash Borer Infestations in Depressional Black Ash Wetlands.” Ecohydrology 10 (7): e1862. https://doi.org/10.1002/eco.1862.

Verdin, Andrew, Balaji Rajagopalan, William Kleiber, and Richard W. Katz. 2015. “Coupled Stochastic Weather Generation Using Spatial and Generalized Linear Models.” Stochastic Environmental Research and Risk Assessment 29 (2): 347–56. https://doi.org/10.1007/s00477-014-0911-6.

Verdin, Andrew, Balaji Rajagopalan, William Kleiber, Guillermo Podestá, and Federico Bert. 2019. “BayGEN: A Bayesian Space-Time Stochastic Weather Generator.” Water Resources Research 55 (4): 2900–2915. https://doi.org/10.1029/2017WR022473.

Watras, C. J., K. A. Morrison, J. L. Rubsam, and I. Buffam. 2017. “Estimates of Evapotranspiration from Contrasting Wisconsin Peatlands Based on Diel Water Table Oscillations.” Ecohydrology 10 (4): e1834. https://doi.org/10.1002/eco.1834.

Watras, C. J., J. S. Read, K. D. Holman, Z. Liu, Y.-Y. Song, A. J. Watras, S. Morgan, and E. H. Stanley. 2014. “Decadal Oscillation of Lakes and Aquifers in the Upper Great Lakes Region of North America: Hydroclimatic Implications: DECADAL WATER LEVEL OSCILLATION.” Geophysical Research Letters 41 (2): 456–62. https://doi.org/10.1002/2013GL058679.

White, Walter Noy. 1932. A Method of Estimating Ground-Water Supplies Based on Discharge by Plants and Evaporation from Soil: Results of Investigations in Escalante Valley, Utah. Vol. 659. US Government Printing Office.

Wilks, D S, and R L Wilby. 1999. “The Weather Generation Game: A Review of Stochastic Weather Models.” Progress in Physical Geography 23 (3): 329–57. https://doi.org/10.1177/030913339902300302.

Wilks, Daniel S. 2012. “Stochastic Weather Generators for Climate-Change Downscaling, Part II: Multivariable and Spatially Coherent Multisite Downscaling: Stochastic Weather Generators for Climate-Change Downscaling.” Wiley Interdisciplinary Reviews: Climate Change 3 (3): 267–78. https://doi.org/10.1002/wcc.167.

Zhu, Jianting, Michael Young, John Healey, Richard Jasoni, and John Osterberg. 2011. “Interference of River Level Changes on Riparian Zone Evapotranspiration Estimates from Diurnal Groundwater Level Fluctuations.” Journal of Hydrology 403 (3-4): 381–89. https://doi.org/10.1016/j.jhydrol.2011.04.016.